Entanglement perturbation theory for the quantum ground states in two dimensions 
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Hilbert space truncation. 

In recent Letters, one of us (SGC) started a novel 
many-body method to calculate partition functions for 
2,3 dimensional macroscopic classical systems and the 
ground states of macroscopic quantum systems in one 
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A simple, general and practically exact method, Entanglement Perturbation Theory (EPT), is 
formulated to calculate the ground states of 2D macroscopic quantum systems with translational 
symmetry. An emphasis will be placed on the applicability of EPT to fermions. We will discuss 
some preliminary evidences which indicate a potential of EPT. 

PACS numbers: 71.10.Fd, 71.27. +a, 75.10.Lp 

To calculate the partition functions and solve the 
Schrodinger equations for macroscopic quantum systems 
are undoubtedly the most fundamental problems in the- 
oretical physics. It would not be an exaggeration to say 
that the history of development of quantum theory is 
primarily that of pursuit of a variety of methods to this 
problem. Methods are roughly categorized into two, a 
variety of mean field theories and otherwise. The essence 
of mean field theories is to make truncations of correla- 
tions |l|. The merit of mean field theories lies in an es- 
sential simplification of calculations. When the systems 
of interest are very complex with many local degrees of 
freedom, such as involving electrons of orbital degeneracy 
and phonons, mean field methods are currently the only 
methods available. The penalty is, nevertheless, concep- 
tually huge: if we don't know the importance of those 
neglected correlations, we are not sure how good the re- 
sults of mean field approximations are. This is precisely 
the starting point and motivation of non-mean field the- 
ories. Here we currently have roughly three methods. 
First are rigorous methods such as Bethe Ansatz 
In spite of its remarkable success over the last several 
decades, a major drawback of rigorous methods is that 
they are mostly limited to 2D classical or ID quantum 
systems of rather special structures. Second are numeri- 
cal simulations such as exact diagonalization and Monte 
Carlo [lOf . With an ever- increasing computer power and 
a steady supply of various ingenious ideas in algorithm, 
they provide a solid starting point often with a superb 
accuracy, but the finite size problem continues to be a 
serious issue. The negative sign problem is yet another 
serious problem in Monte Carlo. Third is the method 
of NRG (numerical renormalization group) , initiated by 



dimension [14|, [15| . Let us call this method as EPT (en- 
tanglement perturbation theory). We here extend the 
EPT to calculate the quantum ground states in two di- 
mensions. Of particular interest is how the EPT handles 
the fermion anticommutation algebra leading to infinite 
range correlations among fermions. Let us call this prob- 
lem simply the sign problem in a wider sense [lif . No ex- 
isting non-mean-ficld methods can handle this sign prob- 
lem satisfactorily. We will show that the 2D fermions 
can be formulated by EPT in a slightly complex but ex- 
ecutable way. We will give some preliminary results for 
the spin 1/2 Heisenberg antiferromagnet and the Hub- 
bard model both on a square lattice, demonstrating a 
significant potential of EPT in concept, no RG and no fi- 
nite size problem, and simplicity and practical exactness 
in numerical implementation. 

Spins and bosons. As an example of spins and bosons, 
let us consider the Heisenberg antiferromagnet (HA) on 
a square lattice. The 2D HA model describes the antifcr- 
romagnetic interaction between near neighbor spins with 
exchange coupling J, 
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where al is the Pauli spin matrix at the site i. To calcu- 
late the ground state of the Schrodinger equation 



(2) 



Wilson [ll|], developed further first by White's DMRG 



(density matrix RG) [HI and most recently by Vidal's 
ER (entanglement renormalization) [13j. However, the 
inability of DMRG for two dimensions is now clear and 
one still needs to see how much that inability of DMRG 
can be improved by ER. Indeed, a deep question remains, 
namely how far one can go with the very idea of the 



we follow the following steps. 

First, instead of ([2]), consider the eigenvalue problem 
for the density matrix 



(3) 



But unlike Monte Carlo and NRG simulations for the 
ground states, our /3 is not the inverse temperature, it 



is a mere parameter here, (3 0, and thus the largest 
eigenstate of the operator 1 — (3H gives the ground state. 
As will become clear below, this simple trick enables EPT 
to efficiently use the power of transfer matrix method. 
Second, we note a decomposition of the density matrix, 
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where "vc" and "he" denote vertical and horizontal 
chains, and the every single chain density matrix is fur- 
ther decomposed into two bond groups: the even group 
connecting the sites (2i, 2% + 1), and the odd group con- 
necting the sites (2i + 1, 2i + 2), where i is the site index 
along the chain, 



-PH cha 



(5) 



Now the local bond density matrix should be further de- 
composed as, 



-l3H h 



1 — f3J(fi ■ dj 
fa ®g a 



(6) 



where and below the repeated indices imply a summa- 
tion, and f a takes four operators, 1, ^/ '/3Ja x , \//3Jo~ y and 
\fl3Jo~ z and g a likewise operators at site j. Thus the ma- 
trix product representation of the even group bonds in 
the density matrix is, ■ • • /„ (g) g a ® ® g@ ® / 7 (g) g 7 • • • , 
and the same expression for the odd group bonds with 
one lattice shifted from the even group case. Putting to- 
gether, we have the matrix representation of the density 
matrix ([5]) for a horizontal chain as, 
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■ Ti p ® rL ® Tlx ® r 2 
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The density matrix for a vertical chain K vc is the same 
as that for a horizontal chain, Kh c , Eq.([7]). Shifting 
one lattice from row to row and column to column, to 
get a checkerboard-type repetition of T 1 ' 2 in horizon- 
tal and vertical directions, we get a 2D extension of 
T 1 ' 2 , namely ri /3 (horizontal)r 2 5 (vertical) -> T^g, and 



r 2 ^(horizontal)r 7(5 (vertical) — > r 2 /37<5 , where we have 
used the same notation T for both ID and 2D. The total 
density matrix in 2D is thus expressed as a bipartite T 1 ' 2 
network as in Fig.l. Note that this is a straightforward 
extension of the ID case to 2D [l5| . 

Third, the ground state wave function can be written 
also as a straightforward extension of the ID case, as in 
Fig.l, where C 1 ' 2 now have 4 legs, left-right-up-down each 
running 1 — p, where p is the entanglement in 2D, in ad- 
dition to the fifth leg representing the 2 local states | f) 
and | 1). This tensor product form of the wave function 
in 2D can be derived by a repeated application of SVD 
(singular value decomposition) just like in ID 



15(. Note 



that the bipartite structure of the total density matrix is 
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FIG. 1: Schematic representation of the density matrix eigen- 
value problem, Eq.(|3} 
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FIG. 2: Schematic representation of the numerator . 



a pretense, the Hamiltonian has a translational symme- 
try. We choose wave function to be of bipartite structure 
to allow for a possible antiferromagnetic phase. We thus 
arrive at a schematic representation, Fig.l, of the density 
matrix eigenvalue problem, Eq.Q. An important com- 
ment here is on a possibility of broken-symmetry ground 
states of different types. For instance, a spin density 
wave-type ground state may be possible for some class of 
systems. The extension of EPT along this line would be 
a very important future problem [l7| . 

Fourth, we consider the variational problem [i = 
^i^vp/vl/vl/ — >. max, where K = e~@ H , by iteration start- 
ing with an input state for ^. The numerator is essen- 
tially a partition function of the 2D network of the ob- 
ject S 1 ' 2 , cf. Fig. 2. The 4 legs of the object S 1,2 , runs 
1 — rp 2 , where r = 4 is for the interaction channel, cf. 
Eq.©. Note that, if S 1 and S 2 are the same, the cal- 
culation is essentially that of the partition function for 
the 2D Ising model [14|. Owing to the macroscopic sys- 
tem size, our job is to calculate the largest eigenvalue 
and eigenvector of a ID infinite object shown in Fig. 3. 
Let us denote this object simply as H. Considering the 
bipartite structure of 3 1,2 and absence of right-left sym- 
metry, the appropriate right eigenvector ipn should have 
a bipartite structure, characterized by £ ' , and likewise 
left eigenvector (p^ characterized by £ 3 ' 4 . Again, we 
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FIG. 3: Schematic representation of a ID infinite object 3 
and its eigenvalue problem. The bottom figure shows a unit 
structure whose repetition along the vertical direction forms 
the quantity ifLEtpR. 



solve this eigenvalue problem by variation. We pursue 
ipL^ l PR/ l PL i PB. ~ > max. A local ingredient of the numer- 
ator is shown in the bottom of Fig. 3. This quantity is 
regarded as a p A r 2 q 2 x p 4 r 2 q 2 matrix A, where q is the 
second entanglement associated with the wave functions 
fL.R- The real nonsymmetric matrix A can be written 
as A = RvL tr where the matrices L, R, and v are made 
up of left eigenvectors, right eigenvectors and eigenval- 
ues of A and tr means the transpose. The eigenvectors 
are normalized as L tr ■ R = 1, and due to this property, 
the summation over the combined entanglement-bond in- 
dices of dimension p 4 r 2 q 2 in the numerator can be done 
N — 1 times, N — > oo in the end, and thus we only keep 
the largest eigenvalue uq and eigenvectors L and R . 



We have ip^^Lpn 



N- 



Lo* r ^4Ro- The denominator 



ifLfR is handled likewise. Let us denote the correspond- 
ing largest eigenvectors as Lo and Ro and eigenvalue as 
p . At this point we realize the "Russian doll" structure, 
namely the above eigenvalue problems and our proce- 
dures to solve them are essentially those in our study of 
the ID model, cf. Fig.l in [15j . The 3 1 ' 2 corresponds to 
r 1,2 in the ID model, £ 1,2 corresponds to f 1,2 in ID, with 
a difference here being the absence of left-right symmetry, 
so that we need the left eigenstate <pl made up of £ 3 ' 4 as 
well as ipn of £ 1 ' 2 . We will not write down rather lengthy 
indexes involved in the resulting eigenvalue equations for 
£ 1,3 and £ 2 ' 4 , but they essentially look like 



(8) 



for i = 1 , 2 where Mi and Ni are square matrices con- 
structed out of 3 1 ' 2 , £ 1-4 , L , Ro, L and Ro- The 
eigenvalue problem for <pl,r is thus reduced to solve these 
two nonlinear, generalized eigenvalue problems by itera- 
tion. Note that the largest eigenvalue kq calculated from 



the two eigenvalue problems should coincide when the 
iteration converges. 

The above procedure is for . We repeat the pro- 

cedure for which is simply to replace T 1 ' 2 in 3 1,2 by 
unity, and the resulting object is denoted as T 1 ' 2 . Corre- 
sponding to fL.R we write Xl,r, corresponding to £ 1-4 
we write f? 1-4 , and we get similar eigenvalue problems 
as Eq.© for ?7 1-4 . Technically, we solve the problem for 
jji— 4 fi rs ^ anc j then use the resulting r? 1-4 as a good in- 
put for £ 1-4 , the reason being that the T 1 ' 2 are close to 
unity with a deviation from it an order of /3 which is very 
small, 10" 6 [H. 

With the converged £ 1-4 and f? 1-4 , we are now ready to 
derive nonlinear, generalized eigenvalue equations for £ > 
resulting from the variational problem, /i = ^fK^f/^f^! — » 
max. Consider first the numerator. Its ID ingredi- 
ent 3 can be written as 3 = ^[k^r where $>l.r are 
the square matrices made up of the left, right eigenvec- 
tors of 3, namely <pl,Ri and k is the diagonal matrix 
made up of the eigenvalues. Note = T r (E N ), and 

use the orthogonal property ^[^r = 1, we then have 
= KQ~ 1 ipi,'Ripji 1 where the suffix "o" indicates 
the largest eigenvalue. On the other hand, the result- 
ing composite is already analyzed above as, tfL^^PR = 
i' v_1 Lo* r ^4Ro- Note that we took in the above varia- 
tion of the last quantity with respect to £ 1_4 . We now 
take variation of the same quantity with respect to C 1 ' 2 . 
The two matrices Mj in Eq.® and Xi below are there- 
fore closely related to each other. The denominator "F^P 
can be analyzed likewise. In the end, the variational 
problem fx = /^V? — > max is reduced to nonlinear, 
generalized eigenvalue problems for C 1 ' 2 . It is written 
schematically as, 



(9) 



for i — 1,2, where Xi and Yi are some symmetric matri- 
ces. We solve this equation for the next C 1 ' 2 until con- 
vergence. 

Finally after the convergence, we can calculate various 
ground state properties. In general, the expectation value 
of the two operators A and B sitting on the adjacent sites 
is calculated as (AB) = 'FAB'F/'F'F. Now 
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h^D(AB)H 



(10) 



where T is obtained from 3 by replacing T 1 ' 2 — > 1, xl.r 
are the left and right eigenvectors of T corresponding to 
the largest eigenvalue erj. Y(AB) is obtained from T by 
inserting AB into a vertical or horizontal bond. Likewise 
the matrix D is obtained from the object at the bottom of 
Fig.3 by replacing £ 1-4 by ry 1 " 4 and 3 1 - 2 by T 1 - 2 . D(AB) 
is then obtained from D by inserting AB into a bipartite 
T 1 ' 2 structure. Ao, Lo and Ro are the largest eigenvalue 
and eigenvectors of D. We thus have 



(AB) 



~D(AB)R /X 



(11) 
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In some preliminary calculations below, we have used 
the parameter (3 = 10 ~ 6 . We have checked (3 = 1CP 4 
with negligible differences. The convergence criterion is 
IIC W " CUII/llCwII < 1 • 10- 8 , and for £ and r,, the 
similar criteria be less than 1-10 7 . When this condition 
is met in the latter case, the relative change in the largest 
eigenvalue often hits 10 -15 , the machine precision. 

Fermions. Let us move on to fcrmions. As a repre- 
sentative 2D fcrmion system, we consider the Hubbard 
model on a square lattice, 



H = —t 



E 

<7,<ij> 



(c 4 W + h.c.) + nttnu (12) 



where t is the transfer integral, U is the onsite Coulomb 
potential and Cj CT , c\ a are the annihilation and creation 
operators for electrons at site i and spin a. We take t 
as the energy unit. To proceed as in spins and bosons, 
we first rewrite the Hamiltonian (TT2"j) as a sum of a local 
bond Hamiltonian, 



bond 



bond 



bond 



with 



Hi = — riifUii - ~(rvr + «u) 



(13) 



(14) 



(15) 



where the chemical potential /i is introduced to control 
the electron number per site. 

Following Eqs.(2-5), Eq.© now becomes 



I 1 - —n^n^ + — (n< t + n 4 )j • [i ->■ j\ 



+I3ty(clc 3a + h.c.) 



f2 Q <8 6 



(16) 

where Vl a takes five operators, 1 — 4r ni T n 4 + ^-{n^ + 
riii), ct^., Cif, ct|_, and C4 and Q likewise operators at 
site j. 

We now come to the point where dimensionality mat- 
ters. We need to specify how we represent the ground 
state. We here order the 2D lattice sites as 

••• 1 2 3 4 5 6 ••• 
•••7 8 9 10 11 12 ••• 
••• 13 14 15 16 17 18 ••• 

where at each lattice site, the four basis states are or- 
dered as |0), I t), I I) and | fT). Let us first consider 
the density matrix associated with the horizontal chains. 
Since the local bond density matrix (fT6l) contains even 
number of creation and annihilation operators, the ma- 
trix representation of the density matrix (1161) reads as 



• • • • • (6 4+ i|(g)(6 l |e- /3ff! "'" d |a. i )(g)|a l+ i) • • '(8)1(8)1 • • • 

where " 1" is unit matrix, and therefore can be written as 
an operator product of local matrices, 



(lk\t 



~ fa,ik ® 9<x,jl 



(17) 



where 
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etc. Note that the —1 in the fa matrix is due to the 
fermion anticommutation algebra. With this new f a and 
g a , Eq.© holds for the fermion Kh c - Note that due 
to our choice of ordering of lattice sites, the fermion 
anticommutation algebra appears only locally for the 
horizontal-chain density matrix. 

The sign problem, a global effect of anticommutation 
algebra, arises for the vertical part. To see this, consider 
the matrix element of a local, vertical-bond density ma- 
trix connecting the sites (i,j) and (i+l,j), where row i is 
one layer above the row i + 1 and column j is one column 
left of column j + 1. The matrix element reads, 

(k+i,j\ ® ® • • • ® (bi, j+ i\ (8 (b i}j \ 

n a (i,j)®e a (i + l,j) 
\ai,j) <8 \a it j + i) (8 • • • (8 |a i+ ij_i) (8 \cn+i,j) 
= f a (i,j) (8 S a (i, j + 1) (8 S a («, j +2) • • • 
• • • ( i + 1 , j - 2 ) (8 (i + 1 , j - 1 ) <8 g Q ( i + 1 , j ) , ( 1 8 ) 

where S a is a unit matrix for the even operator, a = 1 
and 

/ 1 \ 
0-100 
0-10 
V 1 ) 

for the odd operators, a = 2 — 5. 

The operator network of a bipartite F 1 ' 2 with an in- 
finite string of S matrices attached to each and every 
vertical bond, appears hopelessly complex but not im- 
possible. Remember that our parameter /3 — > in the 
density matrix. Thus we only need to retain the j3 linear 
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FIG. 4: Taking variation of £ 1,2 in the row i, indicated by 
an open circle. The unit Hamiltonian is denoted by solid 
lines. The solid circles indicate the fermion T 1 ' 2 while the 
solid rectangles indicate the presence of 5 matrices. 

terms, which brings us back to the original Schrodinger 
equation. From this viewpoint, we can reformulate EPT 
in the following way. Variationally, the Schrodinger equa- 
tion is equivalent to minimize ^2 unit H un u^ />I n F where 
H unit denotes a ID subset of the Hamiltonian as denoted 
by solid lines in Fig. 4. Now due to translational symme- 
try, one can write the numerator as the total number of 
units times an expectation value of H un i t in the row 0. 
Then, a lack of translational symmetry of this quantity 
^Hunit^ in the vertical direction concerning C, 1 ' 2 means 
that we have to variate a vertical string of C , namely 
we need to carry out 

J2^H unit ^) (19) 

i 

where Si means to variate C, 1,2 in the row i, cf. Fig. 4. 
Summing over the entire rows i, the Hamiltonian and 
hence the Schrodinger equation are recovered. Clearly, 
we now have to calculate not only the ground state but 
also the excited states Xl _r, & = 0, 1, 2, • • • of the fermion 
version of the infinite object T and associated eigenval- 
ues. An important note is that, due to the translational 
symmetry of the H u „i t along the horizontal direction, the 
excited states Xl r needed here are only those which have 
the same translational symmetry as the ground state. 
This tremendous simplification is essential in making the 
EPT analysis of 2D fermions practically feasible. In fact, 
recall that T = A#eA^ where are the left and right 
eigenvector xl,r matrices and e the diagonal eigenvalue 
matrix, with orthogonality A^A/j = 1. Inserting this 
expression for T, one can easily see that the main part of 
the calculation in the variation Eq. (fl"9|) is the matrix el- 
ements schematically shown in Fig. 5. We calculate f? 1-4 
as in spins and bosons, but £ 1-4 are not needed here, 
and the variation of Eq. (fl"9]) with respect to C 1 ' 2 leads to 
similar eigenvalue problems as Eq.([9|). 

Preliminary calculations. We have done a preliminary 
calculation on the spin 1/2 HA model, cf. Eq([T|). The 




FIG. 5: Schematic representation of the matrix elements 
which appear in the variation, Eq. (|19l) . 

simplest case p = 1 is soluble by hand. In fact, putting 

C 1 = cos(a)| t) +sin(a)| \) 
C 2 =cos( 7 )|t)+sin( 7 )||), (20) 

we easily find that the energy per bond is \ cos [2 (a — 7)] 
which becomes minimum — j at a — 7 = tt/2, and the 
magnetic moments are icos(2a) and — icos(2a) for the 
sublattices 1 and 2 corresponding to £ 1,2 in the unit of 
the Bohr magneton (1b- Thus the ground state at p = 1 is 
the Neel state with a staggered magnetic moment exactly 
1/2, no quantum reduction. This result was reproduced 
by our EPT to the machine precision. Moving on to the 
case p > 1, we have so far examined the case (p,q)=(3,l) 
with a result; energy per bond is — 0.34J, and the stag- 
gered magnetic moment is 0.39. The energy perfectly 
agrees with the spin wave-Monte Carlo-exact diagonal- 
ization results [Ia420l ]. but the staggered moment 0.39 
is appreciably larger than their result 0.31, rather close 
to the perturbation expansion from the Ising limit to the 
order of J 10 [2lj, which gives 0.385. We certainly need to 
pursue a stable and converging EPT algorithm to reach 
a clear conclusion. 

We have also done a calculation on the 2D Hubbard 
model ignoring the global, infinite string of S matrices 
with an anticipation that this global sign problem should 
become less and less important for larger Coulomb repul- 
sion where electrons are strongly localized. The calcula- 
tion takes much more time, and we have tested only one 
point in the [/-electron number per site plane; U = 8 and 
half- filling. In the case (p, q) = (3, 2), we have found the 
ground state energy per site = —0.473, and a Neel state 
with a staggered magnetic moment 0.45. A remarkable 
precision of EPT can be seen in the sublattice magnetic 
moments, Si and 5*2. We found that Si = —S2 to the 
precision of 10 -4 . In fact, we have chosen the point U = 8 
and half-filling, because a Monte Carlo result was avail- 
able for the energy, — 0.48±0.005 |22j, which is very close 
to our result -0.473. On the other hand, the staggered 
magnetic moment appears to be too large, urging again 
a pursuit of stable and converging EPT algorithm. An 
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interesting point is that, due to the bipartite structure, 
we can check the symmetry of the superconducting order 
parameter, s-wave or d-wave. For instance, the d-wave 
order parameter reads as [23j ]. 

Ad = 71 ^ Cot ° xi ~ c °^ Cx ^ ~ ( c °t c ^i _ c <H<Vr) 
+ (c fC_ x ^ - Coj.C-^) - (cofC-j/J, - CaiC-jrf)], (21) 

where indicates a reference site, ±cc indicates one lat- 
tice to the right or left, and ±y indicates one lattice up 
or down. It should be worth emphasizing that our local 
basis is grand canonical, |0), | f)j I I) an d | an d 
the electron number per site is controlled by the chem- 
ical potential, thus our EPT has an ability to calculate 
the superconducting order parameter. For U — 8 and at 
half-filling with (p,q) = (3,2), we have found no super- 
conductivity, the order parameter |Ad| < 4.10 , as is 
anticipated. 

In conclusion, we have formulated EPT for the quan- 
tum ground states in two dimensions. The key point of 
EPT is that it does not truncate the Hilbert space (no 
RG), nor has the finite size problem. Our preliminary 
calculations on the 2D HA and the Hubbard model show 
that EPT is computationally feasible. A fortunate situa- 
tion for EPT is that, due to the fact that it is based on the 
local basis, the free electron gas is most challenging but of 
course soluble by hand. Moreover, the fermion sign prob- 
lem can still be handled by EPT: An essential difference 
from spins and bosons is that we need to calculate the 
excited states of the ID infinite object T but with an es- 
sential simplification that only the excited states with the 
same translational symmetry as the ground state need to 
be calculated, making the EPT calculation highly exe- 
cutable. We thus recognize the following issues for EPT: 
(1) To establish a stable and converging numerical pro- 
cedure. (2) To calculate the translationally invariant ex- 
cited states of T and treat 2D fermions on a rigorous 
footing. Finally, it is known that Monte Carlo encoun- 
ters a negative sign problem not only in fermions but also 
in spins on a frustrated lattice such as Kagome. EPT is 
free from this difficulty, and the quantum spin liquid issue 
offers an excellent testing ground for EPT. 
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